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Abstract 

In this paper we propose and analyze a new Multiscale Method 
for solving semi-linear elliptic problems with heterogeneous and highly 
variable coefficient functions. For this purpose we construct a general- 
ized finite element basis that spans a low dimensional multiscale space. 
The basis is assembled by performing localized linear fine-scale com- 
putations in small patches that have a diameter of order H log(i7 _1 ) 
where H is the coarse mesh size. Without any assumptions on the 
type of the oscillations in the coefficients, we give a rigorous proof for 
a linear convergence of the H 1 -error with respect to the coarse mesh 
size. To solve the arising equations, we propose an algorithm that is 
based on a damped Newton scheme in the multiscale space. 
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1 Introduction 



This paper is devoted to the numerical approximation of solutions of semi- 
linear elliptic problems with rapidly oscillating and highly varying coefficient 
functions. We are concerned with the following type of equations: 

-V ■ (AVu) + F(u, Vu) = g. 

Here, we prescribe a (zero-) Dirichlet boundary condition for u, A is a highly 
variable diffusion matrix and F is a highly variable nonlinear term that typ- 
ically describes advective and reactive processes. In particular, we have a 
linear term of second order and nonlinear terms of order 1 and 0. A typi- 
cal application is the stationary (Kirchhoff transformed) Richards equation 
that describes the groundwater flow in unsaturated soils (c.f. [IJ[2J[3]). The 
equation reads 

V • (KVu) - V • (Kkr(M(u))g) = f, 

where u is the so-called generalized pressure, K is the hydraulic conductivity 
in the soil, kr the relative permeability depending on the saturation, M is 
some nonlinearity arising from the Kirchhoff transformation and g denotes 
the gravity vector. If we add an infiltration process, the equation receives an 
additional nonlinear reaction term. 

The numerical treatment of such equations is often complicated and ex- 
pensive. Due to the high variability of the coefficient functions, one requires 
extremely fine computational grids that are able to capture all the fine scale 
oscillations. Using standard methods such as Finite Element or Finite Vol- 
ume schemes, this results in systems of equations of enormous size and there- 
fore in a tremendous computational demand that can not be handled in a 
lot of scenarios. 

There is a large variety of so-called multiscale methods that are able to 
overcome this difficulty by decoupling the fine scale computations into local 
parts. This decreases the computational demand without suffering from a 
remarkable loss in accuracy. Examples of multiscale methods are the Het- 
erogeneous Multiscale Method (HMM) by E and Engquist [5] and the Multi- 
scale Finite Element Method (MsFEM) proposed by Hou and Wu [15] . Both 
methods fit into a common framework and are strongly related to numeri- 
cal homogenization (c.f. [HH Efl G31). HMM and MsFEM are typically not 
constructed for a direct approximation of exact solutions but for homoge- 
nized solutions and corresponding correctors instead. This implies that they 
are only able to approximate the exact solution up to a modeling error that 
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depends on the homogenization setting (c.f. [ID])- General proofs of conver- 
gence are therefore hard to achieve. 

We are concerned with a multiscale method that is based on the concept 
of the Variational Multiscale Method (VMM) proposed by Hughes et al. |16j . 
In comparison to HMM and MsFEM, the VMM aims to a direct approxima- 
tion of the exact solution without suffering from a modeling error remainder 
arising from homogenization theory. The key idea of the Variational Mul- 
tiscale Method is to construct a splitting of the original solution space V 
into the direct sum of a low dimensional space for coarse grid approxima- 
tions and high dimensional space for fine scale reconstructions. In this work, 
we consider a modification and extension of this idea that was developed 
in JTHJ [2T] and that was explicitly proposed in [22] • Here, the splitting is 
such that we obtain an accurate but low dimensional space V ms (where we 
are looking for our fine scale approximation instead of an approximation of 
a coarse part) and a high dimensional residual space V { . The construction 
of y ms involves the computation of one fine scale problem in a small patch 
per degree of freedom. Mesh adaptive versions of the VMM with patch size 
control that can be also applied to this multiscale method were achieved in 
[T8| IT9] EHl E3]. The first rigorous proof of convergence was recently obtained 
in (22] for linear diffusion problems. 

In this contribution, we present an efficient way of handling semi-linear 
elliptic multiscale problems in the modified VMM framework, including a 
proof of convergence based on the techniques established in [22]. Even though 
the original problem is nonlinear, the fine scale computations are purely 
linear problems that can be solved in parallel. The main result of this article 
is the optimal convergence of the if 1 -error between exact solution u and its 
multiscale approximation uff. We show that, if the patch size is of order 
Hlog^H^ 1 ), the following error bound holds true: 

\\u-u^\\ HHn) <CH. 

Here, C denotes a generic constant independent of the mesh size of the 
computational grid or the oscillations of A and F. 

The paper is structured as follows. In Section [2] we introduce the setting of 
this paper, including the assumptions on the considered semi-linear problem. 
In Section [3] we present and motivate our method and we state the corre- 
sponding optimal convergence result. The result itself is proved in Section 
HI Finally, to solve the arising nonlinear equations, we propose an algorithm 
that is based on a damped Newton scheme in the multiscale space. This is 
done in Section |5j 
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2 Setting 



Let Q C M. d be a bounded Lipschitz domain with polyhedral boundary, let 
V := Hq(£1) and let A G L°°(Q, denote a matrix valued function with 

uniformly strictly positive eigenvalues. We assume that the space Hq(Q) 
is endowed with the iF-semi norm given by \v \m(n) '■= ||Vu ||i,2(n) (which 
is equivalent to the common iF-norm in Hq(Q)). By (•,•) := (•, -)i 2 (f2) we 
denote the inner product in L 2 (Q) and F : Q xtxM^lisa nonlinear 
measurable function. 

For a given source term g G L 2 (Q) C H^ 1 ^) we are concerned with the 
problem to find u G Hq(Q) (i.e. with a homogeneous Dirichlet boundary 
condition) with 

(AVu,Vv) + (F(;u,Vu),v) = (g,v) (1) 

for all test functions v G Hq(Q). To simplify the notation, we define the 
operator B : H^(Q) ->• if" 1 (ft) by 

(S(u),«;) ff -i tH i := (AVv,Vw) + (F(;v,Vv),w) for v,w E ^(O), 

where (-, -^-i^i denote the dual pairing in if^ft). Here, A is a rapidly 
oscillating, highly heterogeneous diffusion matrix. £) is also allowed 

to be rapidly oscillating, without further assumptions on the type of the 
oscillations. 

We assume that A and F are of the same order of magnitude. This 
prohibits for instance advection dominated processes, which must be treated 
with a different approach and requires a different set of multiscale basis 
functions than the one proposed in this paper. In particular, we will show 
that if the F-term does not dominate the equation, it is sufficient to construct 
a multiscale spaced only based on the oscillations of A. This implies that 
we only need to solve linear problems on the fine scale. For this purpose 
condition (A4) below is to guarantee that A and F have the same size. If 
this is not fulfilled, the proposed method needs modifications with respect to 
the construction of the multiscale basis. Typical examples where condition 
(A4) is not fulfilled are models for transport of solutes in groundwater where 
we deal with extremely large Peclet numbers and a corresponding scaling of 
the advective terms. In this case, the oscillations of F become significant to 
obtain accurate upscaled and homogenized properties (c.f. [TT| IT2]). 

For the subsequent analytical considerations and in order to guarantee a 
unique solution of ([!]), we make the following assumptions: 
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Assumption 1. We assume: 
(Al) A e L°° (tt,R d *£) with 

oo > p := \\A\\ La0 (si) = ess sup sup — — - — . 

xen ceK d \{o} I CI 

and there exists a such that 

(J < a := ess mi mi — — - — , 

x^n C6K d \{0} ICI 

(A2) There exist L\,L% G M>o such that uniformly for almost every x in Q: 

\F(x,S u C)-F(x,6,C)l <^i|Ci-6|, forall(;eR d , 6,6^1, 
|F(x,e,Ci)-i ? (^,e,C 2 )| <L 2 |Ci-C2|, for all CuC2eR d , (GM, 
F(x,0,C) = ^O,£,0) = 0, /or fl«(Gl rf , (el. 

(^45j -B : y — > V is hemicontinuous , i.e. for all Vi,V2,w G Hq(Q) 

s i — y (B(vi + svz),w) H -i H i is a continuous function on [0, 1] 
and B is strongly monotone, i.e. there exist Co > so that for all 

(B(u) - B(v),u-v) H -i iH i > co\u- v\ 2 H1{n) . (2) 

(AJf.) A and F are of similar size, i.e. 0((3)=0(Li)=0(L 2 ). 

Under assumptions (Al)-(A3), the Browder-Minty theorem yields a unique 
solution of problem (PQ). 

Remark 1. Let the C p denote the constant appearing the Poincare-Friedrichs 
inequality for Hq(Q) functions. Observe that (Al)-(A3) imply that the solu- 
tion u G Hl(Q) of (C]j fulfills 

\\F{u, Vn)|| L2(n) < \\F(u, Vu) - F(u,0)\\ L 2 {n) + \\F(u,Vu) - F(0, Vu)||r» ( n) 

< {LxC p + L 2 )|tt| g i(n) < C p LlCp + — \\g\\ L 2(n)- (3) 

Co 
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Note that problem ([TJ also covers equations such as 

-V • (k(u)AVu) = F(«,Vm), 

for a strictly positive and sufficiently regular function k. In this case, the 
equation can be rewritten as 

-V • AVu = F(u,Vu). 

In the remainder of this paper, we use the notation qi < q 2 if q± < Cq 2 
where C > is a constant that only depends on the shape regularity of the 
mesh, but not on the mesh size or the validity of assumption (A4). This 
means that dependencies such as (Li + L 2 )a _1 are always explicitly stated. 
Dependencies on the contrast £ are allowed to be contained in <. 

3 A Rigorous Multiscale Method 

In this section we suggest a Rigorous Multiscale Method (RMM) that is based 
on the concept introduced by Hughes et al. [16] and the specific constructions 
proposed in [181 EI] for linear problems. The required multiscale (MS) basis 
functions are obtained with the strategy established in |22j . 

The main idea of the Variational Multiscale Method is to start from a 
finite element space Vh with a highly resolved computational grid and to 
construct a splitting of this space into the direct sum Vh = V 1 © V f of a 
low dimensional space V 1 and a 'detail space' V f containing all the missing 
oscillations. Then, a basis of V 1 is assembled and we can compute a Galerkin 
approximation ui of « in V 1 . However, the success of this approach strongly 
depends on the choice of V 1 . On the one hand, the costs for assembling 
a basis of V 1 must be kept low. On the other hand, the basis functions 
somehow need to contain information about fine scale features. For instance, 
a standard coarse finite element space is cheap to assemble but will fail to 
yield reliable approximations. On the contrary, the space spanned by high 
resolution finite element approximations yields perfect approximations, but 
is as costly as the original problem that we tried to avoid. Therefore, the 
key is to find an optimal balance between costs and accuracy. In previous 
works (c.f. [HHdHMTH]) ^ ne mu ltiscale basis (MS-basis) of V 1 was constructed 
involving the full multiscale operator B that corresponds with the left hand 
side of the original problem. In a fully linear setting, this can be a reasonable 
choice. However, it gets extremely expensive if B is a nonlinear operator, 
since it leads to numerous nonlinear equations to solve. Furthermore it is not 
clear if the constructed set of basis functions leads to good approximations. 
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One novelty of this work is that we do not involve the full operator B in the 
construction of the MS-basis, but only the linear diffusive part (AV-,V-). 
Even though the oscillations of F are not captured by the MS-basis, we can 
show that we are still able to obtain accurate approximations and to preserve 
the optimal convergence rates. 

3.1 Notation and discretization 

Let 7~h denote a regular triangulation of Q and let H : Q — >■ M >0 denote the 
Tff-piecewise constant mesh size function with H\t = '■— diam(T) for all 
T E Th- Additionally, let Th be a regular triangulation of Q that is supposed 
to be a refinement of Th- We assume that Th is sufficiently small so that all 
fine scale features of B are captured by the mesh. The mesh size h denotes 
the maximum diameter of an element of Th- The corresponding classical 
(conforming) finite element spaces of piecewise polynomials of degree 1 are 
given by 

V H := {v H G C°(0) n Hq(Q) I VT G Th '■ (v h )\t is affine}, 
V h := {v h G C°(Q) n H^Sl) \ VKeT h : (v h )\ K is affine}. 

By J, we denote the dimension of Vh and by Mh = {zj\ 1 < j < J} the 
set of interior vertices of Th- For every vertex Zj G Mh-, let Xj G Vh denote 
the associated nodal basis function (tent function), i.e. Xj G Vh with the 
property Xj(zj) = Sij for all 1 < i,j < J. 

From now on, we denote by Uh G Vh the classical finite element approxi- 
mation of u in the discrete (highly resolved) space Vh, i.e. Uh G Vh solves 

/ AVu h ■ Vv h + F(-, u h , Vu h )v h = / gv h (4) 
Jn Jn 

for all Vh G Vh- We assume that Vh resolves the micro structure, i.e. that 
the error \\u — Uh\\m(n) becomes sufficiently small by falling below a given 
tolerance. For standard finite elements methods the error can be bounded 
by C ■ h s with s > |. However, for regular coefficients, C depends on the 
derivative of A and F with respect to the spatial variable. If A and F are 
rapidly oscillating, the derivative becomes very large and h must be very 
small to equalize the dominance of C. This is only fulfilled, when h resolves 
the micro structure (we refer to [21] and [22] for some qualitative character- 
ization of this so-called resolution condition). We are therefore dealing with 
pre-asymptotic effects for the standard methods. The multiscale method 
that we propose in the subsequent sections will be constructed to converge 



7 



to Ufa with a linear speed in H and with a generic constant that is not af- 
fected by the fine scale oscillations, i.e. in particular we do not have such 
pre-asymptotic effects. 

3.2 Quasi Interpolation 

The key tool in our construction is a linear (quasi-)interpolation operator 
3h '■ Vu — > Vh that is continous and surjective. The kernel of this operator 
is going to be our fine space (or remainder space) V[. In [22J a weighted 
Clement interpolation operator was used. In this work, we do not specify 
the choice. Instead, we state a set of assumptions that must be fulfilled in 
order to derive an optimal convergence result for the constructed multiscale 
method. 

Assumption 2 (Assumptions on the interpolation). We make the following 
assumptions on the interpolation operator 3h '■ Vh — > Vh-' 

(A5) Jh G LiVh^Vn), i-e. Jh is linear and continuous, 

(A6) the restriction ofjjj to Vh is an isomorphism, in particular there holds 
(J H ° ^X^ff) = v H for all v H G V H , 

(A7) there exists a generic constant C^ H , only depending on the shape regu- 
larity ofTn and Tfa such that for all Vh G Vh and for all T G Th there 
holds 

H~ l \\v h - 3 H (vh)\\L2(T) + ||V(% - 3 H {vh))\\v(T) < C 3H \\Vv h \\ L 2 {ulT ) 
with 

Ut ■= \J{k e T H \KnT^0}. 

(A8) there exists a generic constant C^ H , only depending on the shape reg- 
ularity of Th and Th, such that for all vh G Vh there exists Vh G Vh 
with 

3 H (v h ) = v H , Kiffi(n) < Cj H \v H \m(n) and swppv h C suppt^- 

A natural choice for 3h might be the Lagrange interpolation operator 
that fulfills the interpolation estimate in assumption (A7). Furthermore, 
since the Lagrange operator is a projection, the assumptions (A5), (A6) and 
(A8) are automatically fulfilled. 
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Another possibility is to choose 3h as a weighted Clement interpolation 
operator. This construction was proposed in [22J. Given v G Hq(Q), 3hV := 
J2j=i v j^j defines a (weighted) Clement interpolant with nodal values 

v i : = (In vX J dx ) I (In X j dx ) ( 5 ) 

for 1 < j < J (c.f. [SI El E]) and zero in the boundary nodes. Furthermore, 
there exists the desired generic constant Cj H (only depending on the mesh 
regularity parameter and in particular independent of Ht) such that for all 
v G Hq(Q) and for all T G Th there holds 

H~ l \\v - 3 H v\\ L 2 {T) + \\V(v - 3 H v)\\ L 2 {T) < Cj h \\Vv\\ L 2 {u1t) . 

We refer to 0,[6J for a proof of this estimate. This gives us (A7). Assumptions 
(A5) and (A6) are obvious. The validity of (A8) was proved in [22J. 

Note that in certain applications a particular interpolation operator may 
be superior over other choices (cf. Remark 4 in [22] ). 

3.3 Multiscale Splitting and Modified Nodal Basis 

In this section, we construct a splitting of the high resultion finite element 
space Vh i^o a low dimension multiscale space V ms and some high dimen- 
sional remainder space V^. From now on, we let Jh '■ Vh — » Vh denote 
an interpolation operator fulfilling the properties (A5)-(A8). We start with 
defining V£ as the kernel of 3h in Vh- 

Vl := {v h G V h | 3 H v h = 0}. 

V h f represents the features in Vh not captured by Vh- As already mentioned, 
it can be shown that (3h)\v h '■ Vh — > Vh is an isomorphism (see [22]). We 
therefore get 

V h = V H ®Vl where v h = 3 H 1 (3 H (v h )) + v h - 3 H l {3 H {v h ))- (6) 

"-V^ V v ' V v ' 

Here, the property (3h ° ^h^^h) — V H f° r all v h G Vh implies the equation 
3 H (v h - 3 H l (3 H (v h ))) = 3 H (v h ) - (3 H o 3 H 1 )(3 H (v h )) = 0. We still need 
to modify the splitting of Vh, because Vh is an inappropriate space for a 
multiscale approximation. We therefore look for the orthogonal complement 
of Vl in Vh with respect to the inner product (AV-, V-)l2(q). For this purpose, 
we define the orthogonal projection P : Vh — > V h as follows. For a given 
v h G V h , P{v h ) G V h { solves 

(AVP(v h ), Vw f ) = (AVv h , Vw { ) for all w { G V^. 
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Defining the multiscale space V™ h by V™ h := (Vh — P(Vh)), this directly 
leads to an orthogonal splitting: 

V h = Vg% © Vl (7) 

because of 

V h = kern(P) © V l h = (V h -P(V h )) © V h f § (V H -P(V H )) © V< = V™ © V h { . 

Hence, any function Vh G Vh can be decomposed into Vh = v H s +v { with = 
3 H \3 H (v h ))-P(3 H l (3 H (v h ))) and v l = v h - 3 H l (3 H (v h )) + P(3 H 1 (3 H (v h ))). 
Furthermore it holds (AVvh s , Vw f ) = for all w { G V^ f . The space Vff s h is 
a multiscale space of the same dimension as the coarse space Vh- However, 
note that it is only constructed on the basis of the oscillations of A. The 
oscillations of F are not taken into account. We will show that Vff s h still 
yields the desired approximation properties. 

We now introduce a basis of V^ h . The image of the nodal basis function 
Xj G Vh under the fine scale projection P is denoted by <^ = P(Xj) G V h { , 
i.e., 4>j satisfies the corrector problem 

(AV$, Vro) = (AVXj, Vw) for all w G V h { . (8) 

A basis of Vff s h is then given by the modified nodal basis 

{Af := Xj - # | 1 < j < J}. (9) 

As we can see, solving dHJ) involves a fine scale computation on the whole 
domain Q. However, since the right hand side has small support, we are 
able to truncate the computations. As we will see in the next section, the 
correctors show an exponential decay outside of the support if the coarse 
shape function Xj. 

We define a (preliminary) RMM approximation without truncation. 

Definition 2 (RMM approximation without truncation). The Galerkin ap- 
proximation u H s h G Vtf h of the exact solution u of problem (JTJ is defined as 
the solution of 

(AVu^% Vv) + (F(u H %, Vugy , v) = (g, v) for all v G V^ h . (10) 
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3.4 Localization 



So far, in order to construct a suitable multiscale space, we derived a set of 
linear fine scale problems (|S]) that can be solved in parallel. Still, as already 
mentioned in the previous section, these corrector problems are fine scale 
equations formulated on the whole domain Q which makes them almost as 
expensive as the original problem. However, in [22J it was shown that the 
correction <^ decays with exponential speed outside of the support of the 
coarse basis function A^. We specify this feature as follows. Let k G N >0 . 
We define nodal patches Uj^ of k coarse grid layers centered around the node 
zj G Af H by 

:= supp Xj = U {T e T H I Zj e T} , 

r _ , J (11) 

u jtk := U {T eT H \Tn cJ iifc _i ^ 0} for k>2. 

These are the truncated computational domains for the corrector problems 
(jHJ). The fast decay is summarized by the following lemma: 

Lemma 3 (Decay of the local correctors [22]). Let assumptions (Al) and 
(A5)-(A8) be fulfilled. Then, for all nodes Zj G Mh and for all k G N>o there 
holds the following estimate for the correctors 

Remind the definition of '< ' at the end of Section 

This fast decay motivates an approximation of 0^ on the truncated nodal 
patches Uj 7 k- We therefore define localized fine scale spaces by intersecting 
V£ with those functions that vanish outside the patch Uj.k, i.e. 

KW) := {v e Vl | v ]n \^ k = 0} 

for a given node Zj G Nh- The solutions <^ fc G V^(ojj±) of 

(AVfa, Vw) = (AVXj, Vw) for all w G V£(uj jjk ), (12) 

are approximations of <pj from with local support and therefore cheap to 
solve. We define localized multiscale finite element spaces by 

V£J* = spau{A$ := A, - fa | 1 < j < J} C V h . (14) 

We can now define a RMM approximation including truncation: 
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Definition 4 (RMM approximation with truncation). The Galerkin approx- 
imation G of the exact solution u of problem ([T]) is defined as the 
solution of 

(AVu™ k k , Vv) + (F(u™t Vu™ k k ),v) = (g, v) for all v G V^ k . (15) 

Note, that changing F and g does not change the multiscale basis {A™^ | 1 < 
j < J}. Once V^ k is computed, it can be reused for various data functions 
F and g. This makes the new problems cheap to solve. 

3.5 A-priori error estimate 

We are now prepared to state the main result of this contribution, namely 
the optimal convergence of the method for the case that the local patches 
have a diameter of order iJlog(if _1 ): 

Theorem 5. Let u G Hq(Q) denote the exact solution given by problem (T7j) ; 
let Uh G Vh denote the corresponding finite element approximation in the 
Lagrange space with a highly resolved computational grid (i.e. the solution of 
([^]) ) and let VL^' h G be the solution of our proposed multiscale method 

with truncation (i.e. the solution of ( 175)) ). // assumptions (A1)-(A8) are 
fulfilled and if the diameter of the local patches Uj t k fulfills 0(dia,m(ujj^)) > 
H log(||i/ _1 ||ioo(Q)) (1 < j < J), then the following a-priori error estimate 
holds true: 

\\ u ~ u H S ,h\\H^n) < C (\\H\\ Laa{Q ) + \\u - ^Hhi^)) . 

Here, C denotes a generic constant with the property C =0(1). In particular, 
C is independent of the oscillations of A and F and only depends on the 
shape regularity of Th and Th- A suitable choice of m G N with k = m ■ 
log(||// _1 ||l°o(q)) depends on the contrast @-. The larger the contrast, the 
bigger should be m. 

A proof of Theorem[5]is presented in the subsequent section. In particular, 
the result is a conclusion from Theorem [8] which is stated in Section H] below. 
In Theorem |8] we also give details on the generic constant C. We will see 
that it essentially depends on ( Ll + L2 > _ Recall that L\ and L 2 denote the 
Lipschitz constants of F (of. (A2)) and that a is the smallest eigenvalue of 
A. This shows the significance of assumption (A4). For instance, consider 
the scenario of a pollutant being transported by groundwater flow. In this 
case, A describes the hydraulic conductivity which changes its properties on 
a scale of size e. On the other hand, F describes the gravity driven flow 
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that is scaled with the so called Peclet number. However, in the described 
scenario the Peclet number is of order e _1 (c.f. Bourlioux and Majda [I]) 
implying that O(Li) = e _1 . So the generic constant C is of order e -1 . This 
means that we need H < e, i.e. we still need to resolve the micro structure 
with the coarse grid Th producing the same costs as the original problem. 
On the other hand, if H 3> e the estimate stated in Theorem is of no value, 
because the right hand side remains large. 



4 Error Analysis 

This section is devoted to the proof of Theorem In particular, we state a 
detailed version of the result (see Theorem [H] below) , where we specify the 
occurring constants. The proof is splitted into several lemmata. We start 
with an a-priori error estimate for a RMM approximation without truncation: 

Lemma 6. Let Uh G Vh denote the highly resolved finite element approxima- 
tion defined via equation and let u™ h G V™ h denote the RMM approxima- 
tion given by equation [W\l . Under assumptions (Al)-(A3) and (A5)-(A8), 
the following a-priori error estimate holds true: 

\ u h — u'nlmin) 
^ Co(\\Hg\\L 2 (n) + \\H\\ L °o(n)C p 1 P -|blU 2 (n)), 

Co 

where 

+ \\H\\l°°(q.){LiC p +L 2 ) 



C :-- 



Cq ■ a 



Proof. Due to ([7]), we know that there exist u]f h G Vff* h and u h G V£, such 
that 

"ft = + u{. 

We use the Galerkin orthogonality obtained from the equations (jlj) and (TTOj) 
to conclude for all v G Vff s h : 

(AV{u h - «^),Vd) + (F(u h ,Vu h ),v) - (F(u%%Vu™ h ),v) = 0. (16) 
In particular v = u]f h — u^f h G V^ s h is an admissible test function in ([TBI) . 
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Together with 3 H (u h ) = 0, this yields: 

m 

< (AV(u h -u^ h ),V(u h -u^ h )) 

+ (F(u h , Wu h ) - F(u™ h , Vu™ h ),u h - u™ h ) 

® (avk-^.vk-^)) 

= (AVK - «S5i). V«l> + Vtt fc ) - Vtt fc ), 4 - 3h(uD) 

+(F(u%% Vu h ) - F{u^% Vu^ h ), u{ - 3 H {u{)) 

< /3\u h -u^ h \ H i {n) \u[\ H i {n) 

+ \\ H \\L<*>(n)(L 1 \\u h ~ w™fclUa(n) + L 2 \u h - u^ h \ H x {n) )\u{\ H i (n) 

< (/? + ||£r|| L oo ( n)(LiC p + L 2 )) • K - w^l^H") • 



With (AVu^ h , Vu{) = and with 3 H (v { ) = for all v t G 1/ f we get 



= (g,vi-3 H (u { h )) - {F(u h ,Vu h ),u* h -3 H {u { h )) 

13} L C + L 

< (||-H"5lU 2 (n) + ||#lk~(«)C P - -|blU 2 (n)) • 141^(0)- 

Co 

The theorem follows by combing the results. □ 

The next lemma is a conclusion from the previous one: 

Lemma 7. Let G Vh denote the fine scale approximation obtained from 
equation OD and let uf^£ G V^ k denote the solution of problem < T73]) (jfif% 
discrete RMM with truncation). If the assumptions (A1)-(A3) and (A5)- 
(A8) hold true we obtain the following estimate: 

ms.fc I 



Kms./c I 
_ u H,h \HHSi) 

< C 2 \\g\\ L , m \\H\\ L ~ m +C, min JA^(u^ h - v™ h k )\\ LHn) , 



ms.fc ,_ T 7- ms.fc 
»H,h eV H,h 
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where 

' (3 + \\H\\ L o a (Q- ) (LiC p +L 2 ) 



a 



a :=(/? + (LtCp + i^cj- 

h — n i n n ^\C V + L 2 

U 2 ■— C-i + Ui • U p , 

O3 :— . 

Co 

Proof. Let f^f € denote an arbitrary element. Using the Galerkin 

orthogonality obtained from (J3J) and (|T5l) . we start in the same way as in the 
proof of Lemma to get: 

I ms,fc|2 
Co\U h - U Hh 

m k k 

< (AV(u h - u™' h ),V(%- u H ' h )) 

+(F(u h , Vu h ) - F(u™f, Vu™£),u h - u™f) 

@ (Av( Uh - 0,vk - ,,§30 + v(«S5, - «gf )> 

Vu fe ) - , Vu™ h k ), (u h - u%%) + (u%% - v™ h k )) 

< (/3 + [L x C p + L 2 )C p )\u h - u™£\ H i {n) \u h - u™ h \ H i {n) 

+(1 + arifaC, + L 2 )C p )\u h - uSV(n) P*V(u£ fc - ^)IUw 

Dividing by — M^^ln 1 ^) and estimating — M#/J.ffi(fi) with Lemma M 
yields the result. □ 

Combining the results of Lemma [3] and Lemma [7] yields the main result 
of this contribution: 

Theorem 8. Let Uh £ Vh be solution of and let u™ ,k G V^ k be the 
solution of 173]) . If the assumptions (Al)-(A3) and (A5)-(A8) hold true and 
if the number of layers k fulfills k > logdl.ff' -1 !!^^)), then it holds 



(n), 



where 



C :— C 2 + C p — C3 

Co 



and u«£/i C2 and C3 as in Lemma^ 
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Proof. We define w™£ G Vg% by 



j j 

ms,fc 



W H,h 



3=1 i=i 



where u™ h (zj), j = 1, 2, . . . , J, are the coefficients in the basis representation 
of 

u lih fr° m Definition [2J Hence, 



min J|^V(^-<?)ll! 2( n) 



<\\A^(u^-w-f)\\l 2{n) 



For details on the last step, we refer to Lemma 18 in [22]. Due to the Galerkin 
orthogonality for the corrector problems we get 

ll^ 1/2 V(0? - ^ k )\\h {n) < ||^v#||£^ Wfc _ . 

We refer to the first part of the proof of Lemma 8 in [22] for details. The 
application of Lemma [31 (JED and some inverse inequality yield 

ll^ 1/2 V(0j-0j )fc )||| 2(n) < e- 2( ^ )1/2fe ||^V0j||i 2(n) 

< e- 2 ^ 1/2k \\A 1 / 2 VX'j\\l, iu) 

< /3e- 2 ^ 1/2fc i/- 2 ||Aj||i 2(n) 

By choosing k = m ■ log(||i7 _1 ||z°°(n)) with m e N, we can achieve an ar- 
bitrary fast polynomial convergence of this term in H (this will also cancel 
the k d term). However, we bound this by a linear convergence since this is 
fastest rate that we can obtain for the whole error. Finally, the combination 
of this estimate and (HH) plus £/ =1 ^(%) 2 ||Aj||i 2(n) < p H ^||f 2(n) < 

\\^ u H%\\h { n) < C p c o 2 \\9\\l^n) y ields the assertion. □ 

5 The Multiscale Newton scheme 

In this section we discuss a solution algorithm for handling the nonlinear 
multiscale problem ( TT5|) . For this purpose, we consider a damped Newton's 
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method in the multiscale space V^ k . Remind the considered problem: we 
are looking for u £ Hq with 

(B(u),v) H -i jH i = (g,v) for all v £ H%(Sl), 

where we introduced the notation 

(B(v),w) H -^ H i := (AVv,Vw) + (F(-,v,Vv),w). 

Here, B : Hq(Q) — > is a hemicontinuous and strongly monotone 

operator due to assumption (A3). As already mentioned, under these as- 
sumptions, the Browder-Minty theorem yields a unique solution of the above 
problem. However, we will need an additional assumption on F to guarantee 
that the Newton scheme converges: 

Assumption 3. Let DF(x, •, •) denote the Jacobian matrix of F(x, •, •). 

(A9) We assume that there exists some constant Lo > so that for almost 
every x in Q and for all (£, Q £ R x M. d 

\DF(x, 6, Ci) - DF(x, 6, C 2 )| < L D \(^, Ci) - (6, C 2 )|, 

i.e. F(x,-,-) £ W 2 '°°(R x R d ). 

For clarity of the presentation we will leave out several indices within this 
section. In particular, we make use of the following notation: 

Definition 9. For simplicity, we define 

V ms ■= Vffl with basis Af := \% = Xj - <P), k for 1 < j < J. 
Furthermore, we denote u ms := u^'£ . Additionally, let 

d 1 F(x,CX)--=d ( F(x^,C) and d 2 F(x, C) := d c F(x, £, C). 

We now describe the Newton strategy in detail. The fully discrete multi- 
scale problem to solve reads: 

find u ms £ V ms : (AVu ms , VXf s ) + (F{-, u ms , Vu ms ),Xf s ) - (g, Xf s ) = 

for all 1 < j < J. Again, using Browder-Minty, u ms exists and is unique. 
Accordingly, we get the following well posed algebraic version of the problem: 

find a £ R J : G(a) = 
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and where G : IR J — > M J is given by 

(G(a)\ (18) 
j j j 

:= ^(AVAf , VAD + (/••(••^n,Af.^n,VAf).An - (g^™). 
j=i j=i j=i 

We have the relation w ms = X^f=i ^'A™ s . Before we can apply the Newton 
method to (|18|) . we need to ensure that the iterations of the scheme are well 
defined. The following lemma ensures this: 

Lemma 10. Let (X, || • ||x) denote a Hilbert space with dual space X' . Let 
furthermore B : X — >■ X' be a hemicontinuous and strongly monotone oper- 
ator on X, i.e. there exists Co > so that 

(B(v) — B(w),v — w)x > Cq\\v — w\\ 2 x for all v,w G X and 
s i — y {B{u + sv),w)x 

is a continuous function on [0, 1] for all u,v,w G X . Let Xn denote a finite 
dimensional subspace with basis {tpi, ...,iJjn} and let b : W N — > Vn define the 
linear bijection with b(a) := *Yld =l ctiipi- If G(a) := b~ l (A(b(a))) , then the 
Jacobi matrix DG(a) G ~R nxn has only positive eigenvalues. 

Proof. Let B' denote the Frechet the derivative of B, given by 

r>// W \ i- B{u + sv) - B{u) 

B (u)(v) = hm — — for u,v G X. 

s->0 S 

This and the strong monotonicity yield: 

,„„ w s v , (B(u + sv) - B(u))(v) 

u s— >0 S 

= lxm—(B(u+sv) - B(u))(u + sv -u) (19) 

s->0 S 2 

> lim 4tC ||sw|| 2 = c ||u|| 2 . 

s-i>0 s 

Next, observe that b induces an inner product on M. N by («i, a.2)b '■= (b(ai), b(a.2))x- 
Let a := b~ l {u) then we get 

B(u)m = hrn ^ s 

v (b o IrWEjUte + Aytyj) - (b o i- 1 )^ a^)) 
= lim 

s^O S 

G(aj + sej) — G(ctj 



b ( i im 1 ~-Zi± 

. s-tO S 



b(D a G{a)ej) 
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Using this, we get for arbitrary £ 6 M. N and v% := fr(£) : 

JV 

(D a G(a)e,e)6 = X>&(A*G(a)ei,e;) 6 

iV 

= ^60(^aG(a)e,),6(e,))x 
= X)^(B'(u)(^),^)x 

in 

= (B r { U ){vt),v € ) x > coHWl = coU\\l 
Since all norms in ~R N are equivalent we have the desired result. □ 

Now, we can apply the damped Newton method for solving the nonlinear 
algebraic equation G(a) = 0. If D a G denotes the Jacobian matrix of G, we 
get the following iteration scheme: 

c^+i) : = a (n) + Aa {n \ 

where Aa'"' solves 

D a G(a (n) )Aa {n) = -G(a (n) ). (20) 
Here, D a G is given by: 

J J 

D ai (G(a)\ := (AWXT S , VAH + (dtF^ajXf^a^X^Xr, Af> 

i=i j=i 

j j 
+ (d 2 F(-, £ a,Af , £ OjVAJ 18 ) • VAr, Aj"">. 

3=1 3=1 

Lemma [TU1 ensures that equation f[2"Uj) has a unique solution Aa'™', i.e. that 
the Newton iteration is well posed. Since G G C 1 (M Ar ) has a nonsingu- 
lar Jacobian matrix D a G (due to Lemma fTUi) and since we have Lipschitz- 
continuity of D a G (due to Assumption [3]), we have the that Newton scheme 
converges quadratically as long as the starting value is close enough to the 
exact solution (c.f. [8J). However, this means that we can only guarantee 
local convergence of the method. In order to ensure global convergence, we 
can use a simple damping strategy, i.e. we are looking for a damping param- 
eter C G (0, 1] so that a( n+1 > := + (Aa^ with the property G(a^ n+1 '>) < 
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G(a^). Under the same assumptions (i.e. (Al)-(A3) and (A9)), we get that 
the damped Newton scheme converges, i.e. there exists a (damping) interval 
(0, Co] and for any ( G (0, Co] we get the above 'damping property'. This is 
an easy observation if we look at the function h(Q := \G(a^ + (Aa^)\ 2 
which fulfills h(0) > and h'(0) = -2G(a^) ■ G(qW) < 0. In summary we 
have globally linear convergence of the method (using damping) and locally 
(i.e. in an environment of the solution) even quadratic convergence using the 
classical Newton scheme without damping. 

With these considerations, we can state the full algorithm below. Recall 
that Nh denotes the set of interior vertices of Th and for Zj G A/#, Xj G Vh 
denotes the corresponding nodal basis function. 

Note that in the presented algorithm, each iteration starts with the damp- 
ing parameter X n = 1 and we do not use damping parameters from previous 
iterations. The advantage is that we automatically get quadratic conver- 
gence of the Newton scheme as soon as we leave the region where damping 
is required. Therefore, damping is only used when really necessary. 

Proposition 11. We use the notation stated in Definition^ Let u G Hq(Q) 
denote the solution of let Uh G Vh denote the solution of O) and let u ms G 
V ms denote the solution of $15\)). Furthermore, we let u ms '^ := u H s ' h k '^ 
define the n 'th iterate from the damped Newton Variational Multiscale Method 
stated in the algorithm. Under assumptions (A1)-(A9), the Newton step Kffl^l 
is well posed, yields an unique solution and u ms '^ converges at least linearly 
to u ms . If furthermore k > log(||i^ _1 ||z,°°(n)), the following a-priori error 
estimates hold true: 



where C is a generic constant with the property C =0(1) (see Theorem^ 
andlBfor details) and 



Here, we have L n (H) < 1. 

jj u ms,(n-i) £ S sufficiently close to u ms , we even get quadratic convergence 
of the Newton scheme, i. e. : 



u - u ms \\ H i {n) < C (\\H\\ L «,(n) + \\u - u h \\ H i {n) ) 



u 



ms ms,(n) 



H i (n) < L n (H)\\u ms - u 



ms,(ri— 1) 



u 



ms ms,(n) 



m(n) < L n (H)\\u ms - u 



ms,(n— 1) ||2 



with 



\\(D a G) 



-i 




L n (H) < 



(R N ) 



L 
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Algorithm: dampedNewtonRMM( abstol, reltol, a^°\ k ) 



In parallel foreach zj G Mh do 

compute 4> h ^ k G Vl{u jtk ) with 

(AV<j)^ k , Vw) = (AVXj, Vw) for all w G V£(u jtk ). 
end 

Set V£5 fc := {A, - | 1 < j < J}. Set A£ = A, - # >Jfc . 

Set qW := a(°). Set <^' (n) := ^ ,nf^. Set 

J J J 

(0(a)), := X^VA^VAE) + <^£«;A£>£«;VA£) 
j'=i j'=i i=i 

Set tol : = |G(a (0) )| 2 • reltol + abstol. 

while |G(a^)| 2 > tol do 
' Set^^^E^r^A-l. 

Define the entries of the stiffness matrix by: 

:= {AVX^ VXtl) + (d 1 F(.,u^ n \Vu^)X^ X^) 
+(9 2 F(,^ ) V<' w )-VA-,A->. 

Define the entries of the right hand side by: 

p( n ) ._ I n \ms\ / /lY7„, ms . fc :( n ) Y7\ms\ / 77./ ms,fc,(n) y ms,fc,(n)\ x 

Find (Aa) (n+1) G M J , with 

M (n) (Aa) (n+1) = F (n) . 

Set A n := 1. Set := a» + X n Aa^ n \ 

while > G(a^) do 

I Set A n := |A n . Set «( n+1 ) := a» + A n Aa< n >. 
end 

Set «W := a> +1 ). Set to/ := |G(a( n ))| 2 • re/toZ + abstol. 
end 

Q„ f „ ms,fe,(n) ^ J ™( n ) \ms 



ms\ 
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where L denotes the Lip schitz- constant of D a G. As indicated, L n (H) typi- 
cally depends on the mesh size. However, in some cases of semi-linear prob- 
lems, it is possible to bound L n (H) independent of the triangulation (c.f. 
\nj ). In particular , if F(x,u,Vu) = F(x,u) (i.e. no dependency on Vu) we 
get that L n (H) = L n independent of the underlying mesh. The proof can be 
obtained analogously to the proof of Proposition 4-1 in WH - The proof fails 
for general F(x, u, Vu) . 

Remark 12. Note that the proposed method only requires the computation 
of the multiscale basis {A™ s | 1 < j < J} once at the beginning. For each iter- 
ation step of the damped Newton scheme, Kffl^l is a cheap (low dimensional) 
linear problem that can reuse the initially computed multiscale basis. If the 
multiscale basis was computed using the nonlinear term F , local corrector 
problems would have to be solved for each Newton step newly, making the 
whole procedure expensive. 

References 

[1] H. W. Alt and S. Luckhaus. Quasilinear elliptic-parabolic differential 
equations. Math. Z., 183(3):311-341, 1983. 

[2] H. Berninger. Non-overlapping domain decomposition for the Richards 
equation via superposition operators. Lect. Notes Comput. Sci. Eng., 
70: Springer, Berlin, 169-176, 2009. 

[3] H. Berninger, R. Kornhuber and O. Sander. On nonlinear Dirichlet- 
Neumann algorithms for jumping nonlinearities. Domain decomposition 
methods in science and engineering XVI, Lect. Notes Comput. Sci. Eng. 
(55), Springer, Berlin: 489-496 (2007). 

[4] A. Bourlioux and A. J. Majda, An elementary model for the valida- 
tion of flamelet approximations in non-premixed turbulent combustion. 
Combust. Theory Model. 4(2): 189-210, 2000. 

[5] C. Carstensen, Quasi-interpolation and a posteriori error analysis in 
finite element methods. M2 AN Math. Model. Numer. Anal, 33 (6): 1187- 
1202, 1999. 

[6] C. Carstensen and R. Verfiirth. Edge residuals dominate a posteriori 
error estimates for low order finite element methods. SIAM J. Numer. 
Anal, 36(5):1571-1587 (electronic), 1999. 



22 



[7] Ph. Clement. Approximation by finite element functions using lo- 
cal regularization. Rev. Frangaise Automat. Informat. Recherche 
Operationnelle Ser., RAIRO Analyse Numerique, 9(R-2):77-84, 1975. 

[8] J. E. Dennis Jr. and R. B. Schnabel Numerical Methods for Uncon- 
strained Optimization and Nonlinear Equations. SIAM Classics in Ap- 
plied Mathematics (1996). 

[9] W. E and B. Engquist. The heterogeneous multiscale methods. Com- 
mun. Math. Sci, 1(1):87-132, 2003. 

[10] A. Gloria. An analytical framework for the numerical homogenization of 
monotone elliptic operators and quasiconvex energies. SIAM Multiscale 
Model. Simul, 5(3):996-1043 (electronic), 2006. 

[11] P. Henning and M. Ohlberger. The heterogeneous multiscale finite el- 
ement method for advection-diffusion problems with rapidly oscillating 
coefficients and large expected drift. Netw. Heterog. Media, 5 (4): 711- 
744, 2010. 

[12] P. Henning and M. Ohlberger. A Note on Homogenization of Advection- 
Diffusion Problems with Large Expected Drift. Z. Anal. Anwend., 
30(3):319-339, 2011. 

[13] P. Henning and M. Ohlberger. A-posteriori error estimation for a het- 
erogeneous multiscale method for monotone operators and beyond a pe- 
riodic setting. FB 10 , Universitdt Miinster, Preprint 04/11 - N, 2011. 

[14] P. Henning. Convergence of MsFEM approximations for elliptic, non- 
periodic homogenization problems. Netw. Heterog. Media, 7(3):503-524, 
2012. 

[15] T. Hou and X.-H. Wu. A multiscale finite element method for elliptic 
problems in composite materials and porous media. J. Comput. Phys., 
134(1):169-189, 1997. 

[16] T.J.R. Hughes, G.R. Feijoo, L. Mazzei and J.-B. Quincy. The vari- 
ational multiscale method - a paradigm for computational mechanics. 
Comput. Methods Appl. Mech. Engrg., 166: 3-24, 1998. 

[17] J. Karatson. Characterizing Mesh Independent Quadratic Convergence 
of Newton's Method for a Class of Elliptic Problems. J. Math. Anal., 
44(3):12791303, 2012. 



23 



[18] M. G. Larson and A. Malqvist. Adaptive variational multiscale methods 
based on a posteriori error estimation: energy norm estimates for elliptic 
problems. Comput. Methods Appl. Mech. Engrg., 196(21-24) :2313-2324, 
2007. 

[19] M. G. Larson and A. Malqvist. An adaptive variational multiscale 
method for convection-diffusion problems. Comm. Numer. Methods En- 
grg., 25(l):65-79, 2009. 

[20] M. G. Larson and A. Malqvist. A mixed adaptive variational multi- 
scale method with applications in oil reservoir simulation. Math. Models 
Methods Appl. Sci., 19(7):1017-1042, 2009. 

[21] A. Malqvist. Multiscale methods for elliptic problems. Multiscale Model. 
Simul, 9(3):1064-1086, 2011. 

[22] A. Malqvist and D. Peterseim. Localization of Elliptic Multiscale Prob- 
lems. larXiv:1110.0692V 3 

[23] J. M. Nordbotten. Adaptive variational multiscale methods for multi- 
phase flow in porous media. SI AM Multiscale Model. Simul., 7(3): 1455- 
1473, 2008. 

[24] D. Peterseim and S.A. Sauter. Finite Elements for Elliptic Problems 
with Highly Varying, Non-Periodic Diffusion Matrix. SIAM Multiscale 
Model. Simul, 10(3):665-695, 2012. 

[25] D. Peterseim. Robustness of Finite Element Simulations in Densely 
Packed Random Particle Composites. Netw. Heterog. Media, 7(1) : 113- 
126, 2012. 



24 



